/* (c) Copyright 1993 by Joseph Felsenstein.
   Written by Joseph Felsenstein, Akiko Fuseki, Sean Lamont, and Andrew Keeffe.
   Permission is granted to copy and use this program provided no fee is
   charged for it and provided that this copyright notice is not removed. */

#include <boost/shared_ptr.hpp>
#include <boost/shared_array.hpp>
#include <math.h>
#include <iostream>
#include <sstream>
#include "KITSCH_nearest_neighbor.h"
#include "tree_to_newhampshire_factory.h"
#include "tree_to_paralogy_factory.h"
#define DEBUG_LEVEL 5
#include "debug_func.h"

using std::vector;
using std::string;
using std::cerr;
using std::cout;
using std::ostringstream;
using boost::shared_array;
using paralogy_tree::t_tree;
using paralogy_tree::t_node;

using boost::shared_ptr;
using boost::shared_array;



namespace kitsch_tree
{

	
//	for (unsigned i = 1; i < cnt_leaves; ++i)
//		for (unsigned j = 0; j < i; ++j)
//	for (unsigned i = 0; i < cnt_leaves - 1 ; ++i)
//		for (unsigned j = i + 1; j < cnt_leaves; ++j)
//	two forms are identical
	
	


//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	struct t_kitsch_node
		
//
enum eCHILD_NUM{CHILD_LEFT, CHILD_RITE, CHILD__ALL};
#define epsilon         0.000001   /* a very small but not too small number */

struct t_kitsch_node
{		  /* describes a tip species or an ancestor */
	t_kitsch_node(unsigned count) : index(0), tip(false), d(count, 0.0),
						w(count, 1.0), leaf_length(0.0), sametime(false), weight(0.0), 
						processed(false),
						parent(0), child_num(CHILD_LEFT)
	{
		clear();
	}
	void				clear ()				{	parent = child[0]  = child[1] = 0;}
	unsigned index;
	bool tip;
	vector<double> d;
	vector<double> w;		/* distances and weights                  */
	double leaf_length;
	bool sametime;			/* bookkeeps scrunched nodes              */
	double weight;			/* weight of t_kitsch_node used by scrunch         */
	bool processed;			/* used by evaluate                       */
	t_kitsch_node *parent;	// should be equivalent to back for tips and main nodes
	t_kitsch_node *child[2];	// should be equivalent to next->back for main nodes
							// should be equivalent to next->next->back for main nodes
	eCHILD_NUM		child_num;
	void				CHECK_LINKS() const;
	t_kitsch_node*		sibling() const;
	void				set_parent(t_kitsch_node* new_parent, eCHILD_NUM num);
	t_kitsch_node& 		operator = (const t_kitsch_node& other);
	bool 				operator!= (const t_kitsch_node& other) const	
															{return !(*this == other);}
	bool 				operator== (const t_kitsch_node& other) const;
};

struct t_kitsch_tree
{
	t_kitsch_tree(	vector<shared_array<double> >&  d,
					vector<shared_array<double> >&  w);
				
	vector<shared_array<double> >&	distance;
	vector<shared_array<double> >&	weights;
	unsigned 									cnt_seq;
	unsigned 									cnt_nodes;
	unsigned									cnt_trees_examined;
	vector<kitsch_tree::t_kitsch_node>			shared_mem;
};

//

//	struct t_kitsch_node
		
//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888



}



//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
// globals
namespace
{

	vector<kitsch_tree::t_kitsch_node* > 				treenode, besttree;
	kitsch_tree::t_kitsch_node*						root		= 0;
	kitsch_tree::t_kitsch_node*						best_root	= 0;
	
	unsigned	njumble				= 5;
	double		power				= 2.0;
	unsigned	cnt_trees_examined	= 0;
	double		like				= 0.0;
	double		bestyet				= 0.0;
	inline	double square(const double d){return d * d;}
}
//
//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888




namespace kitsch_tree
{


	
	
//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	struct t_kitsch_node
		
//


//________________________________________________________________________________________
		
//	CHECK_LINKS
	
//________________________________________________________________________________________
void t_kitsch_node::CHECK_LINKS() const
{
#ifndef NDEBUG
	static unsigned indent = 0;
	if (!indent)
		debug_cout(5, "CHECK_LINKS ENTRY");
	debug_cout(5, string(indent, ' ') << index);
	indent += 2;
	if (parent)
		assert(parent->child[child_num] == this);
	if (!tip)
	{
		if (child[CHILD_LEFT])
		{
			assert(child[CHILD_LEFT]->parent == this);
			assert(child[CHILD_LEFT]->child_num == CHILD_LEFT);
			child[CHILD_LEFT]->CHECK_LINKS();
		}
		if (child[CHILD_RITE])
		{
			assert(child[CHILD_RITE]->parent == this);
			assert(child[CHILD_RITE]->child_num == CHILD_RITE);
			child[CHILD_RITE]->CHECK_LINKS();
		}
	}
	indent -= 2;
	if (!indent)
	{
		for (unsigned i = 0; i < treenode.size(); ++i)
		{
			t_kitsch_node* a = treenode[i];
			assert(a);
			if (a->parent && a->parent->child[a->child_num] != a)
			{
				cout << debug_prefix << a->index  << " failed " << (a->child_num == CHILD_LEFT ? 
							 "left\n" : "right") << std::endl;
				cout << debug_prefix  << a->parent->index  << " parent"  << std::endl;
				if (a->parent->child[a->child_num])
					cout << debug_prefix  << a->parent->child[a->child_num]->index << " child" << std::endl;
				assert(a->parent->child[a->child_num] == a);
			}
			if (a->child[CHILD_LEFT])
			{
				if (a->child[CHILD_LEFT]->parent != a)
				{
					cout << debug_prefix  << a->index  << " failed" << std::endl;
					cout << debug_prefix << a->child[CHILD_LEFT]->index << std::endl;
					if (a->child[CHILD_LEFT]->parent)
						cout << debug_prefix << a->child[CHILD_LEFT]->index << std::endl;
					assert(a->child[CHILD_LEFT]->parent == a);
				}
				if (a->child[CHILD_LEFT]->child_num != CHILD_LEFT)
				{
					cout << debug_prefix << a->index  << " failed" << std::endl;
					cout << debug_prefix << a->child[CHILD_LEFT]->index << std::endl;
					cout << debug_prefix << (a->child[CHILD_LEFT]->child_num == CHILD_LEFT ? 
							 "left\n" : "right") << "should be CHILD_LEFT" << std::endl;
					assert(a->child[CHILD_LEFT]->child_num == CHILD_LEFT);
				}
			}
			if (a->child[CHILD_RITE])
			{
				if (a->child[CHILD_RITE]->parent != a)
				{
					cout << debug_prefix << a->index  << " failed" << std::endl;
					cout << debug_prefix << a->child[CHILD_RITE]->index << std::endl;
					if (a->child[CHILD_RITE]->parent)
						cout << debug_prefix << a->child[CHILD_RITE]->index << std::endl;
					assert(a->child[CHILD_RITE]->parent == a);
				}
				if (a->child[CHILD_RITE]->child_num != CHILD_RITE)
				{
					cout << debug_prefix << a->index  << " failed" << std::endl;
					cout << debug_prefix << a->child[CHILD_RITE]->index << std::endl;
					cout << debug_prefix << (a->child[CHILD_RITE]->child_num == CHILD_RITE ? 
							 "RITE\n" : "right") << "should be CHILD_RITE" << std::endl;
					assert(a->child[CHILD_RITE]->child_num == CHILD_RITE);
				}
			}
	
		}
		debug_cout(5, "CHECK_LINKS EXIT");
	}
#endif
}


//________________________________________________________________________________________
		
//	sibling
	
//________________________________________________________________________________________
t_kitsch_node* t_kitsch_node::sibling() const
{
	assert (parent);
	return parent->child[(child_num + 1) % CHILD__ALL];
}
//________________________________________________________________________________________
		
//	set_parent
	
//________________________________________________________________________________________
void t_kitsch_node::set_parent(t_kitsch_node* new_parent, eCHILD_NUM num)
{
	debug_func_cout(5);
	// old parent's child should be null
	if (parent)
	{
		assert(parent->child[child_num] == this);
		parent->child[child_num] = 0;
	}
	
	parent = new_parent;
	
	if (new_parent)
	{
		// new parent's old child should be orphaned
		if (new_parent->child[num])
			new_parent->child[num]->parent = 0;
		new_parent->child[num] = this;
	}
	child_num = num;
}


//________________________________________________________________________________________
		
//	assignment operator =
	
//________________________________________________________________________________________
t_kitsch_node& t_kitsch_node::operator = (const t_kitsch_node& other)
{
	child_num	= other.child_num;
	d			= other.d;
	w			= other.w;
	leaf_length	= other.leaf_length;
	sametime	= other.sametime;
	weight		= other.weight;
	processed	= other.processed;
	index		= other.index;
	tip			= other.tip;
	return *this;
}

//________________________________________________________________________________________
		
//	operator==
	
//________________________________________________________________________________________
bool t_kitsch_node::operator== (const t_kitsch_node& other) const
{
	return (
			d			== other.d				&&
			w			== other.w				&&
			leaf_length	== other.leaf_length	&&
			sametime	== other.sametime		&&
			weight		== other.weight			&&
			processed	== other.processed		&&
			index		== other.index);
}
//

//	struct t_kitsch_node
		
//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888















//________________________________________________________________________________________

//	initialize_tree

//________________________________________________________________________________________
void initialize_tree(	vector<kitsch_tree::t_kitsch_node>& shared_mem,
						unsigned& next_mem,
						vector<kitsch_tree::t_kitsch_node*>& tree,
						unsigned cnt_seq, unsigned cnt_nodes)
{
	debug_func_cout(1);
	
	// all pointers set to zero
	tree.assign(cnt_nodes, 0);
	
	// all pointers set to allocated memory and set index of nodes
	// indices are immutable
	for (unsigned i = 0; i < cnt_nodes; ++i)
	{
		assert(next_mem < shared_mem.size());
		tree[i] = &shared_mem[next_mem++];
		tree[i]->index = i;
	}

	// mark leaves
	for (unsigned i = 0; i < cnt_seq; ++i)
		tree[i]->tip = true;



}

//________________________________________________________________________________________

//	tree_alloc_mem

//________________________________________________________________________________________
void tree_alloc_mem(vector<kitsch_tree::t_kitsch_node>& shared_mem,
					unsigned cnt_seq, unsigned cnt_nodes)
{
	debug_func_cout(1);
	
	// allocate memory for nodes
	shared_mem.assign((cnt_seq + (cnt_seq - 1)) * 2, t_kitsch_node(cnt_nodes));
	unsigned next_mem = 0;
	
	// assign memory for both current tree and the "best tree" with the lowest
	//		mean square distance
	initialize_tree(shared_mem, next_mem, treenode, cnt_seq, cnt_nodes);
	initialize_tree(shared_mem, next_mem, besttree, cnt_seq, cnt_nodes);
	assert(next_mem == shared_mem.size());
}  /* tree_alloc_mem */


//________________________________________________________________________________________

//	assigndata

//________________________________________________________________________________________
void assigndata(std::vector<boost::shared_array<double> >& matrix,
					unsigned cnt_seq, unsigned cnt_nodes)
{
	debug_func_cout(1);
	
		
	// assign data to nodes distances and weights
	for (unsigned i = 0; i < cnt_seq - 1 ; ++i)
		for (unsigned j = i + 1; j < cnt_seq; ++j)
		{
			if (matrix[i][j] < epsilon)
				treenode[i]->d[j] = 
				treenode[j]->d[i] = epsilon;
			else
				treenode[i]->d[j] = 
				treenode[j]->d[i] = matrix[i][j];
			if (matrix[i][j] > 5000.0 || matrix[i][j] < 0.0)
				treenode[i]->w[j] = 
				treenode[j]->w[i] = 0.0;
			else
				treenode[i]->w[j] = 
				treenode[j]->w[i] =	1.0 / exp(2.0 * log(treenode[i]->d[j]));
		}
	
	
	
}



//________________________________________________________________________________________

//	add

//	inserts the nodes new_parent and its left descendant, newtip to the tree.  
//	at becomes new_parent's right descendant

//________________________________________________________________________________________
void add(t_kitsch_node* at, t_kitsch_node* new_tip, t_kitsch_node* new_parent)
{
#ifndef NDEBUG
	debug_func_cout(1);
	assert(at);
	assert(new_tip);
	assert(new_parent);
	if (at->parent)
		debug_cout(5, at->parent->index << " to new parent");
	debug_cout(5, new_parent->index << " add top");
	debug_cout(5, new_tip->index << " left");
	debug_cout(5, at->index << " right");
	root->CHECK_LINKS();
#endif
	
	
	
	// make new_parent child of at's parent in its place
	new_parent->set_parent(at->parent, at->child_num);
	
	// now link up at and new_tip as children of new_parent
	new_tip->set_parent(new_parent, CHILD_LEFT);					
	at->set_parent(new_parent, CHILD_RITE);
	
	
	// reroot if necessary
	if (root == at)
		root = new_parent;
	root->parent = 0;							




#ifndef NDEBUG
	root->CHECK_LINKS();
#endif

	
}  /* add */




//________________________________________________________________________________________

//	remove

//	removes nodes item and its ancestor, parent, from the tree.
//		the new descendant of parent's parent (in the place of parent) 
//		will be item's sibling
//		returns pointers to the deleted nodes, item and parent

//________________________________________________________________________________________
void remove(t_kitsch_node *& item, t_kitsch_node*& parent)
{
#ifndef NDEBUG
	debug_func_cout(1);
	root->CHECK_LINKS();
	assert(item);
	debug_cout(5, item->index << " remove item");
	debug_cout(5, item->parent->index << " and remove parent");
	if (item->parent->parent)
		debug_cout(5, item->parent->parent->index << " from grandparent");
	if (item->sibling())
		debug_cout(5, item->sibling()->index << " sibling will be in the place of item");
#endif
	

	// forget it if has no parent: e.g. root
	parent = item->parent;						
	if (!parent)						
		return;                                 
	
	t_kitsch_node* sibling = item->sibling();
	
	
	
	// reroot if necessary to the other child of root
	if (root == parent)
		root = sibling;
	if (sibling)
		sibling->set_parent(parent->parent, parent->child_num);
	
	// break up links between item and parent as well
	item->parent		=0;
	parent->child[0]	=
	parent->child[1]	=0;

	
	
#ifndef NDEBUG
	assert(!parent->parent);
	root->CHECK_LINKS();
#endif
	
}  /* remove */




//________________________________________________________________________________________

//	scrunchtraverse

//________________________________________________________________________________________
void scrunchtraverse(t_kitsch_node *u, t_kitsch_node*& closest, double *tmax)
{
	debug_func_cout(4);
#ifndef NDEBUG
	cout << u->index << " scrunchtraverse" << std::endl;
	assert(tmax);
	assert(u);
#endif
	
	//cout << "scrunch traverse tmax=" << *tmax << std::endl;
	/* traverse to find closest t_kitsch_node to the current one */
	if (!u->sametime)
	{
		if (u->leaf_length > *tmax)
		{
			closest = u;
			*tmax = u->leaf_length;
		}
		return;
	}
	if (u->parent)
		u->leaf_length = u->parent->leaf_length;
	if (!u->tip)
	{
		scrunchtraverse(u->child[CHILD_LEFT], closest, tmax);
		scrunchtraverse(u->child[CHILD_RITE], closest, tmax);
	}
}  /* scrunchtraverse */




//________________________________________________________________________________________

//	combine

//________________________________________________________________________________________
void combine(t_kitsch_node& a, t_kitsch_node& b)
{
#ifndef NDEBUG
	debug_func_cout(4);
#endif

	
	/* put t_kitsch_node b into the set having the same time as a */
	if (a.weight + b.weight <= 0.0)
		a.leaf_length = 0.0;
	else
		a.leaf_length = (a.leaf_length * a.weight + b.leaf_length * b.weight) / 
						(a.weight + b.weight);
	a.weight += b.weight;
	b.sametime = true;
}  /* combine */




//________________________________________________________________________________________

//	scrunch

//________________________________________________________________________________________
void scrunch(t_kitsch_node& s)
{
	debug_func_cout(3);
	/* see if nodes can be combined to prevent negative lengths */
#ifndef NDEBUG
	s.CHECK_LINKS();
	root->CHECK_LINKS();	
#endif
	t_kitsch_node *closest	= 0;
	do
	{
		double tmax	= -1.0;
		if (!s.tip)
		{
			scrunchtraverse(s.child[CHILD_LEFT], closest, &tmax);
			scrunchtraverse(s.child[CHILD_RITE], closest, &tmax);
		}
		// found = (tmax > s.leaf_length);
		if (tmax <= s.leaf_length)
			break;
		assert(closest);
		combine(s, *closest);
	} while (1);
}  /* scrunch */




//________________________________________________________________________________________

//	secondtraverse

//________________________________________________________________________________________
void secondtraverse(t_kitsch_node* a, 
					t_kitsch_node* q, 
					t_kitsch_node* u, 
					t_kitsch_node* v, 
					unsigned i, unsigned j, unsigned k, double *sum)
{
	debug_func_cout(4);
#ifndef NDEBUG
	cout << a->index << std::endl;
	assert(a);
	assert(sum);
	assert(q);
	assert(u);
	assert(v);
#endif
	/* recalculate distances, add to sum */

	if (!(a->processed || a->tip))
	{
		secondtraverse(a->child[CHILD_LEFT], q,u,v,i,j,k,sum);
		secondtraverse(a->child[CHILD_RITE], q,u,v,i,j,k,sum);
		return;
	}
	if (!(a != q && a->processed))
		return;
	unsigned l = a->index;
	double wil = u->w[l];
	double wjl = v->w[l];
	double wkl = wil + wjl;
	double wli = a->w[i];
	double wlj = a->w[j];
	double wlk = wli + wlj;
	q->w[l] = wkl;
	a->w[k] = wlk;
	if (wkl <= 0.0)
		q->d[l] = 0.0;
	else
		q->d[l]	= (wil * u->d[l] + wjl * v->d[l]) / wkl;
	if (wlk <= 0.0)
		a->d[k] = 0.0;
	else
		a->d[k]	= (wli * a->d[i] + wlj * a->d[j]) / wlk;
	if (wkl > 0.0)
	{
		double TEMP = u->d[l] - v->d[l];
		(*sum) += wil * wjl / wkl * (TEMP * TEMP);
	}
	if (wlk > 0.0)
	{
		double TEMP = a->d[i] - a->d[j];
		(*sum) += wli * wlj / wlk * (TEMP * TEMP);
	}
}  /* secondtraverse */




//________________________________________________________________________________________

//	firstraverse

//________________________________________________________________________________________
void firstraverse(t_kitsch_node *q_, t_kitsch_node* r, double *sum)
{
	debug_func_cout(3);
  /* firsttraverse                              */
   /* go through tree calculating branch lengths */
	/* Local variables for firstraverse: */
	//cout << "f sum=" << *sum << std::endl;
	if (q_ == 0)
		return;
	
	t_kitsch_node *q = q_;
	q->sametime = false;
	if (!q->tip)
	{
#ifndef NDEBUG
		cout << q->index << " [1]\n";
		assert(q->child[CHILD_LEFT]);
		cout << q->child[CHILD_LEFT]->index << " [1]\n";
		assert(q->child[CHILD_RITE]);
		cout << q->child[CHILD_RITE]->index << " [1]\n";
#endif
		firstraverse(q->child[CHILD_LEFT], r, sum);
		firstraverse(q->child[CHILD_RITE], r, sum);
	}
	assert(q);
	q->processed = true;
	if (q->tip)
		return;
	t_kitsch_node * u = q->child[CHILD_LEFT];
	t_kitsch_node * v = q->child[CHILD_RITE];
	assert(u);
	assert(v);
	unsigned i = u->index;
	unsigned j = v->index;
	unsigned k = q->index;
#ifndef NDEBUG
	cout << i << " [2]" << std::endl;
	cout << j << " [2]" << std::endl;
	cout << k << " [2]" <<std::endl;
	assert(i < v->w.size());
	assert(i < v->d.size());
	assert(j < u->w.size());
	assert(j < u->d.size());
#endif
	if (u->w[j] + v->w[i] <= 0.0)
		q->leaf_length = 0.0;
	else
		q->leaf_length = (u->w[j] * u->d[j] +
				v->w[i] * v->d[i]) /
			   (2.0 * (u->w[j] + v->w[i]));
	q->weight = u->weight + v->weight + u->w[j] + v->w[i];
//	if (!negallowed)
	assert(q);
   	scrunch(*q);
#ifndef NDEBUG
	r->CHECK_LINKS();
#endif
	secondtraverse(r,q,u,v,i,j,k,sum);
}  /* firstraverse */




//________________________________________________________________________________________

//	sumtraverse

//________________________________________________________________________________________
double sumtraverse(t_kitsch_node *q)
{
	debug_func_cout(3);
	
	/* traverse to finish computation of sum of squares */
	assert(!q->tip);
	double sum = 0.0;
	if (!q->child[CHILD_LEFT]->tip)
		sum += sumtraverse(q->child[CHILD_LEFT]);
	if (!q->child[CHILD_RITE]->tip)
		sum += sumtraverse(q->child[CHILD_RITE]);
	t_kitsch_node * u = q->child[CHILD_LEFT];
	t_kitsch_node * v = q->child[CHILD_RITE];
	unsigned i = u->index;
	unsigned j = v->index;
	double TEMP = u->d[j] - 2.0 * q->leaf_length;
	double TEMP1 = v->d[i] - 2.0 * q->leaf_length;
	sum += u->w[j] * (TEMP * TEMP) + v->w[i] * (TEMP1 * TEMP1);
	return sum;
	//cout << "t sum=" << *sum << std::endl;
}  /* sumtraverse */




//________________________________________________________________________________________

//	evaluate

//		fill in times and evaluate sum of squares for tree

//________________________________________________________________________________________
void evaluate(t_kitsch_node& r)
{
	debug_func_cout(1);
	

	for (unsigned i = 0; i < treenode.size(); ++i)
		treenode[i]->processed = treenode[i]->tip;
	
	double sum = 0.0;
	firstraverse(&r, &r, &sum);
	if (!r.tip)
		sum += sumtraverse(&r);
	++cnt_trees_examined;
	like = -sum;
}  /* evaluate */





//________________________________________________________________________________________

//	tryadd

//		temporarily adds one fork and one tip to the tree.
//		if the location where they are added yields greater
//		"likelihood" than other locations tested up to that
//		time, then keeps that location as there

//________________________________________________________________________________________
void tryadd(t_kitsch_node& p,  t_kitsch_node*& item,  t_kitsch_node*& nufork, t_kitsch_node*& there)
{
#ifndef NDEBUG
	debug_func_cout(1);
#endif
	
	add(&p, item, nufork);
	assert(root);
	evaluate(*root);
	if (like > bestyet)
	{
		bestyet = like;
		there = &p;
	}
	remove(item, nufork);



#ifndef NDEBUG
	root->CHECK_LINKS();
#endif
}  /* tryadd */




//________________________________________________________________________________________

//	addpreorder

//		traverses a binary tree, calling tryadd
//		at a t_kitsch_node before calling tryadd at its descendants

//________________________________________________________________________________________
t_kitsch_node* addpreorder(t_kitsch_node& p, t_kitsch_node* item, t_kitsch_node* nufork, t_kitsch_node* there)
{
#ifndef NDEBUG
	debug_func_cout(1);
//	assert(p);
//	if (p == 0)
//		return there;
#endif
	
	tryadd(p, item, nufork, there);
	if (!p.tip)
	{
		assert(p.child[CHILD_LEFT]);
		assert(p.child[CHILD_RITE]);
		there = addpreorder(*p.child[CHILD_LEFT], item, nufork, there);
		there = addpreorder(*p.child[CHILD_RITE], item, nufork, there);
	}

	return there;
}  /* addpreorder */




//________________________________________________________________________________________

//	tryrearr

//		evaluates one rearrangement of the tree.
//		if the new tree has greater "likelihood" than the old
//		one sets success := TRUE and keeps the new tree.
//		otherwise, restores the old tree

//________________________________________________________________________________________
void tryrearr(t_kitsch_node *p,  t_kitsch_node*& r, bool *success)
{
#ifndef NDEBUG
	debug_func_cout(1);
#endif	

	t_kitsch_node * parent		= p->parent;
	if (!parent)
		return;
	
	t_kitsch_node * grandparent	= parent->parent;
	if (!grandparent)
		return;
	
	t_kitsch_node * sibling		= p->sibling();
	
	
	double oldlike = like;
	
	remove(p, parent);
	add(grandparent, p, parent);
	if (r->parent)
		r = r->parent;
	assert(r);
	evaluate(*r);
	if (like > oldlike)
	{
		bestyet = like;
		*success = true;
		return;
	}
	remove(p, parent);
	add(sibling, p, parent);
	if (r->parent)
		r = r->parent;
	like = oldlike;

#ifndef NDEBUG
	root->CHECK_LINKS();
#endif
}  /* tryrearr */




//________________________________________________________________________________________

//	repreorder

//________________________________________________________________________________________
void repreorder(t_kitsch_node *p,  t_kitsch_node*& r, bool *success)
{
	debug_func_cout(1);
	/* traverses a binary tree, calling PROCEDURE tryrearr
	   at a t_kitsch_node before calling tryrearr at its descendants */
	if (p == 0)
		return;
	tryrearr(p, r, success);
	if (!p->tip)
	{
#ifndef NDEBUG
		debug_cout(5, p->index);
		debug_cout(5, p->child[CHILD_LEFT]->index << "left");
		debug_cout(5, p->child[CHILD_RITE]->index << "right");
#endif
		repreorder(p->child[CHILD_LEFT], r,success);
		repreorder(p->child[CHILD_RITE], r,success);
	}
}  /* repreorder */




//________________________________________________________________________________________

//	rearrange

//________________________________________________________________________________________
void rearrange(t_kitsch_node *& r_)
{
	debug_func_cout(1);
	/* traverses the tree (preorder), finding any local
	   rearrangement which decreases the number of steps.
	   if traversal succeeds in increasing the tree's
	   "likelihood", PROCEDURE rearrange runs traversal again */
/* Local variables for rearrange: */
	t_kitsch_node *r = r_;
	bool success = false;
	do
	{
		success = false;
		repreorder(r, r, &success);
	} while (success);
}  /* rearrange */







//________________________________________________________________________________________

//	dtraverse

//________________________________________________________________________________________
/*void dtraverse(t_kitsch_node *q)
{
	// print table of lengths etc.
	unsigned i;

	if (!q->tip)
		dtraverse(q->next->back);
	if (q->back != 0)
	{
		fprintf(outfile, "%4ld  ", q->back->index - cnt_seq);
		if (q->index <= cnt_seq)
		{
			for (i = 0; i < namelength; i++)
				putc(names[q->index][i], outfile);
		}
		else
			fprintf(outfile, "%4ld      ", q->index - cnt_seq);
		fprintf(outfile, "%13.5f", treenode[q->back->index]->leaf_length - q->leaf_length);
		fprintf(outfile, "%15.5f\n", root->leaf_length - q->leaf_length);
	}
	if (!q->tip)
		dtraverse(q->next->next->back);
}  dtraverse */




//________________________________________________________________________________________

//	describe

//________________________________________________________________________________________
void describe(double& sum_sqrs, double& avg_percent_stddev, unsigned cnt_seq)
{
	debug_func_cout(1);
	/* prints table of lengths, times, sum of squares, etc. */
	sum_sqrs = -like;
	double totalnum = 0.0;
	for (unsigned i = 0; i < cnt_seq; ++i)
	{
		for (unsigned j = 0; j < cnt_seq; ++j)
			if (i != j && treenode[i]->d[j] > 0.0)
				totalnum += treenode[i]->w[j] * treenode[i]->d[j] * treenode[i]->d[j];
	}
	totalnum -= 2;
	avg_percent_stddev = 100 * sqrt(-(like / totalnum));
}





//________________________________________________________________________________________

//	copy_kitsch_tree

//________________________________________________________________________________________
void copy_kitsch_tree(vector<t_kitsch_node*> from, vector<t_kitsch_node*> to,
					unsigned cnt_seq, unsigned cnt_nodes)
{
	debug_func_cout(1);
	/* make a copy of a tree */
	for (unsigned i = 0; i < cnt_nodes; ++i)
	{
		*to[i] = *from[i];
		if (from[i]->parent)
		{
			// restore link to parent
			to[i]->parent = to[from[i]->parent->index];
			to[i]->parent->child[to[i]->child_num] = to[i];
		}
		else to[i]->parent	= 0;
	}
}




//________________________________________________________________________________________

//	maketree

//________________________________________________________________________________________
void maketree(std::vector<boost::shared_array<double> >& matrix, double& best_likelihood,
				const vector<unsigned>& order_of_entry,
				unsigned cnt_seq, unsigned cnt_nodes,
				t_progress_indicator& dots)
{
	debug_func_cout(1);
	for (unsigned i = 0; i < cnt_nodes; ++i)
		treenode[i]->clear();
	/* constructs a binary tree from the pointers in treenode.
	   adds each t_kitsch_node at location which yields highest "likelihood"
	   then rearranges the tree for greatest "likelihood" */
	
	// first two nodes
	root = treenode[order_of_entry[0]];
	add(treenode[order_of_entry[0]], treenode[order_of_entry[1]], treenode[cnt_seq]);
	
	// everything else
	for (unsigned i = 2; i < cnt_seq; ++i)
	{
		++dots;
		bestyet = -10000000.0;
		t_kitsch_node* item = treenode[order_of_entry[i]];
		t_kitsch_node * nufork = treenode[cnt_seq - 1 + i];
		t_kitsch_node* there = addpreorder(*root, item, nufork, 0);
		add(there, item, nufork);
		like = bestyet;
		rearrange(root);
		assert(root);
		
		// evaluate without bumping up the number of trees examined
		evaluate(*root);
		--cnt_trees_examined;
	}
	
	
	double bestlike = bestyet;
	double gotlike;
	do
	{
		gotlike = bestlike;
		t_kitsch_node * nufork = treenode[cnt_nodes - 1];
		for (unsigned j = 0; j < cnt_nodes; ++j)
		{
			bestyet = -32000.0;
			t_kitsch_node* item = treenode[j];
			if (item != root)
			{
				remove(item, nufork);
				t_kitsch_node* there = addpreorder(*root, item, nufork, root);
				add(there, item, nufork);
			}
		}
	} while (bestlike > gotlike);
	
	debug_cout(4, "bestlike " << bestlike);
	if (bestlike > best_likelihood)
	{
		copy_kitsch_tree(treenode, besttree, cnt_seq, cnt_nodes);
		best_root = besttree[root->index];
		best_likelihood = bestlike;
#ifndef NDEBUG
		root->CHECK_LINKS();
		best_root->CHECK_LINKS();
#endif
	}
	
	
}  /* maketree */





//________________________________________________________________________________________

//	treeout

//________________________________________________________________________________________
template<typename TREEMAKER> void kitsch_treeout(	TREEMAKER& maker,
													t_kitsch_node *p, 
													const vector<string>& names)
{
	debug_func_cout(1);
	
	assert(p);
	if (p->tip)
	{
		maker.add_leaf(names[p->index]);
	}
	else
	{
		maker.add_internal_node();
		kitsch_treeout(maker, p->child[CHILD_LEFT], names);
		maker.add_sibling();
		kitsch_treeout(maker, p->child[CHILD_RITE], names);
		maker.reparent_to_grandparent();
	}
	if (p == root)
	{
		maker.finish();
	}
	else
	{
		maker.set_branch_len(p->parent->leaf_length - p->leaf_length);
	}
}
													
string kitsch_tree_to_str(const vector<string>& names)
{
	newhampshire_tree_maker make_newhampshire_tree;
	kitsch_treeout(make_newhampshire_tree, root, names);
	return make_newhampshire_tree.str;
}


//----------------------------------------------------------------------------------------

//	kitsch_make_tree

//----------------------------------------------------------------------------------------
void kitsch_make_tree(	paralogy_tree::t_tree& tree,
						const std::vector<std::string>& names,
						std::vector<boost::shared_array<double> >& matrix,
						double& sum_sqrs, double& avg_percent_stddev, 
						unsigned& count_trees_examined,
						t_progress_indicator& dots)
{
	debug_func_cout(1);
	
	// number of genes and nodes
	unsigned cnt_seq = matrix.size();
	unsigned cnt_nodes = cnt_seq * 2 - 1;
	
	// allocate memory for tree
	vector<kitsch_tree::t_kitsch_node> shared_mem;
	tree_alloc_mem(shared_mem, cnt_seq, cnt_nodes);
	assigndata(matrix, cnt_seq, cnt_nodes);
	
	// which order leaves are added
	vector<unsigned> order_of_entry;
	for (unsigned i = 0; i < cnt_seq; ++i)
		order_of_entry.push_back(i);
	
	
	cnt_trees_examined = 0;
	double best_likelihood = -1e50;
	for (unsigned i = 0; i < njumble; ++i)
	{
		maketree(matrix, best_likelihood, order_of_entry, cnt_seq, cnt_nodes, dots);
		for (unsigned i = 0; i < cnt_seq; ++i)
			std::swap(order_of_entry[i], order_of_entry[rand() % cnt_seq]);
	}
	
	// use best tree
	copy_kitsch_tree(besttree, treenode, cnt_seq, cnt_nodes);
	assert(best_root);
	root = treenode[best_root->index];
#ifndef NDEBUG
	root->CHECK_LINKS();
	best_root->CHECK_LINKS();
#endif
	evaluate(*root);
	
	describe(sum_sqrs, avg_percent_stddev, cnt_seq);
	count_trees_examined = cnt_trees_examined;


	paralogy_tree_maker make_paralogy_tree(tree);
	kitsch_treeout(make_paralogy_tree, root, names);
	dots.finish();
}

//________________________________________________________________________________________

//		kitsch_make_tree

//			wrapper

//________________________________________________________________________________________
void kitsch_make_tree(	paralogy_tree::t_tree& tree,
						const std::vector<std::string>& names,
						std::vector<boost::shared_array<double> >& matrix,
						double& sum_sqrs, double& avg_percent_stddev, 
						unsigned& cnt_trees_examined)
{
	debug_func_cout(1);
	ostringstream ignore_strm;
	t_progress_indicator dots(ignore_strm, 10000);
	kitsch_make_tree(tree, names, matrix, sum_sqrs, avg_percent_stddev, 
						cnt_trees_examined, dots);
}






}

